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Abstract. We present a numerical algorithm to solve the Boltzmann equation for the electron distribution 
function in magnetic multilayer heterostructures with non-collinear magnetizations. The solution is based 
on a scattering matrix formalism for layers that are translationally invariant in plane so that properties 
only vary perpendicular to the planes. Physical quantities like spin density, spin current, and spin-transfer 
torque are calculated directly from the distribution function. We illustrate our solution method with a 
systematic study of the spin-transfer torque in a spin valve as a function of its geometry. The results 
agree with a hybrid circuit theory developed by Slonczewski for geometries typical of those measured 
experimentally. 
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1 Introduction 



Magnetic multilayer structures have attracted a great deal 
of experimental and theoretical attention. One motivation 
for these studies is the potential application of such struc- 
tures for data storage and spin dependent transistors. One 
special effect in a magnetic multilayer that can induce 
magnetic reversal or magnetic dynamics is called spin- 
transfer. In the ten years since its theoretical prediction, 
|H2j spin-transfer has been studied extensively both ex- 
perimentally [3I4I5I6I7] and theoretically. |8I9I10I11I12H5] 
One fundamental issue is how to reliably calculate the 
spin-transfer torque in spin valve systems. To solve this 
problem, different approaches have been developed, in- 
cluding the Boltzmann equation [14115] , microscopic quan- 
tum mechanics [16], drift-diffusion theory [17I18I7I19T20] . 
and circuit theory [21|22j . 

Each approach has its own advantages and disadvan- 
tages. Simple theories, for example circuit theory or the 
drift-diffusion approach, treat the transport in terms of 
densities and current densities and do not track individ- 
ual electrons. Such methods have the advantage that they 
can give analytic results in some limits and the disad- 
vantage that they leave out some essential physics. One of 
the major approximations of these theories is that they ig- 
nore the differences between electrons propagating in dif- 
ferent directions. Fully quantum mechanical calculations 
track all of the electrons and all of the coherent scattering 
processes like coherent multiple scattering between layers. 
However, such calculations are quite time consuming and 
the coherent multiple scattering between layers that is in- 



cluded in such an approach does not seem to play a role 
in experimental results. 

The semiclassical Boltzmann equation is a useful com- 
promise between these extremes. In such an approach, 
the scattering is treated semiclassically. For transport in 
collinear magnetic systems, it has been used in two ways. 
Superlattices that have mean free paths longer than layer 
thicknesses can be treated as artificial bulk materials [23124] . 
This approach retains the coherent multiple scattering be- 
tween the interfaces. The other approach is to solve the 
Boltzmann equation within layers and join the solutions 
through boundary conditions at the interface. Here, we de- 
scribe a generalization of the latter approach to treat non- 
collinear magnetizations. In this approach, the Boltzmann 
equation tracks individual electrons through the distribu- 
tion function, but ignores the coherent multiple scattering. 
It is easier to treat defect scattering in such an approach 
than it is in a fully coherent calculation. The advantage of 
the Boltzmann calculation is that it is simply computable 
and includes the essential physics. One disadvantage is 
that it cannot give analytical results. The neglect of co- 
herent multiple scattering between interfaces is both an 
advantage and a disadvantage. The greater simplicity that 
results in the calculation is an advantage for the metallic 
devices that have been measured to date because the co- 
herent scattering does not appear to play a role. However, 
this neglect would be a disadvantage in devices in which 
such effects were important. 

Slonczewski has developed a hybrid approach com- 
bining aspects of circuit theory with a simplified Boltz- 
mann equation to give an analytic expression for the spin- 
transfer torque [25] . This approach gives much more ac- 
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curate results than the drift-diffusion method for typical 
device geometries. However, it breaks down when layer 
thicknesses become longer than are typical. 

In this paper, we describe a numerical algorithm that 
solves the Boltzmann equation for the electron distribu- 
tion function in a magnetic multilayer system using a scat- 
tering matrix formalism. The spin-transfer properties are 
calculated from the resulting distribution function. We 
compare our Boltzmann results to the results from more 
approximate methods. These comparisons show Slonczewski 
hybrid theory is highly accurate for typical experimental 
structures, but drift-diffusion shows systematic deviations. 
We have reported results of calculations using the methods 
described in the present paper in earlier papers [14115] . 

The paper is organized as follows: Section[2]gives a gen- 
eralized spin-dependent Boltzmann equation appropriate 
for ferromagnetic systems; Section [3] describes in detail 
the algorithm that solves the Boltzmann equation in a 
spin valve magnetic multilayer system; Section U shows 
the results of the Boltzmann method developed in Sec- 
tion [3] applied to a spin valve. Section [5] gives a summary. 
Readers who arc not interested in the formalism can skip 
directly to Section [U 



2 Generalized matrix Boltzmann equation 

In non-magnetic materials, the spin independent Boltz- 
mann equation is 



where the integration over k' is now a two-dimensional 
integral restricted to the Fermi surface. The scattering 
rate Pk,k' has been rescaled. A delta function has been 
factored out of each term. In the second term on the left 
hand side, this delta function comes from dfo/de. 

For spin dependent magnetic materials, one needs a 
spin dependent distribution function as well as a spin de- 
pendent Boltzmann equation. When there is a natural 
quantization axis, i.e., all electron spins are either parallel 
f{a =|) or anti-parallel (a =J.) to the axis, an electron dis- 
tribution function is separated into the distribution func- 
tions for spin- up and spin-down electrons: f*(r, k) and 
/l(r,k) 



dff(r,k) dg(r,k) 
v k 5 ehrVk- 



dr 



9ei 



= / dk'P k ,k'[s(r,k')-<?(r,k)], 



where g(r, k) is an electron distribution function that de- 
pends on the spatial coordinate r and electron wave- vector 
k. Vk is the electron velocity, ek its energy, E is the elec- 
tric field, and we have ignored the magnetic field. Pk,k' is 
the probability of electron scattering from state k' to state 
k. The Boltzmann equation is valid only when the bulk 
properties vary slowly. It cannot be used for abrupt inter- 
faces or boundaries; later we describe how to use boundary 
conditions to relate solutions of the Boltzmann equation 
in different regions across the interfaces. 

Transport in metals is dominated by the electrons near 
the Fermi energy. The occupancy of states far from the 
Fermi energy does not change and those states do not 
contribute to the transport. This suggests the use of the 
linearized Boltzmann equation in which the distribution 
function is assumed to have the form 



5 (r,k) = /o(k) + /(r,k)<S(£ F -e k ), 



(2) 



where /o is the equilibrium distribution function. The 
delta function restricts the wave vector to the Fermi sur- 
face, |k| = fcp for free electrons. With this approximation 
for the distribution function, the Boltzmann equation be- 
comes 



d /(r,k) 
dr 



3E-v k = f dk^P k , k ,[/(r,k')-/(r,k)], 

JFS 

(3) 



v * pF v 

Vk • — 5 eE • Vk 

or 



/ dk'p k v[r(k')-r(k)] 

JFS 

+ 1 rfk'p^ k ,[r'(k')-r(k)], 

Jfs 



(4) 



The r dependence of the distribution function has been 
suppressed, and a =|, J, and a ^ a'. Compared with 
Eq. © , Eq. ^ has an additional spin flip scattering term 
on the right hand side because the distribution function in 
Eq. ([3]) includes both spin types, while the one in Eq. (j4| 
is for only one spin type. Similar to the definition of Pk,k' 
in Eq. ([3]), P k k / and P k f k / are the probabilities of electron 
scattering from state k' to state k without and with spin 
flip. We assume that P k f k , is the same for spin flip in both 
directions, up to down or down to up. 

In a non-magnet, there is no natural quantization axis, 
so it is convenient to use the axis of neighboring fer- 
romagnetic layers if the magnetizations of those layers 
are collinear. In this case, the scattering Pk,k' is spin- 
independent, so the substitutions f° = ^(/^ + /I) and 
f z = h(f ~ / ) gi ye the pah of equations 



Vk • 



/°(k) 
dr 



eE-Vk - / dk'P k ,k'[/°(k')-/ (k)] 

JFS 



Vk 



+ / dk'ps f k ,[/°(k')-/°(k)], 

JFS 

= f dk'P k . k '[r(k')-/ 2 (k)] 
ar JFS 

- 1 dk'p£ k ,[r(k')+/ z (k)], 

JFS 

(5) 



In the first equation, spin flip scattering acts as another 
form of non-flip scattering as far as the number accumu- 
lation is concerned. In the second equation, the electric 
field plays no role because it does not couple to the spin 
accumulation. In the systems of interest here, the magne- 
tizations are not collinear, so the spin axis in the distri- 
bution function can vary with both r and k. However, the 
generalization is straightforward. The second equation in 
Eq. is replicated for each of the other directions in spin 
space, with z — > x, y. The generalization can be derived 
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Fig. 1. Schematic view of a perpendicular spin valve structure. 



from a matrix form of the Boltzmann equation in terms 
of the matrix distribution function 



bulk layer in the spin valve; (|3.4p construct the interface 
scattering matrix for each interface in the spin valve; (|3.5[) 
connect the bulk and interface scattering matrices into a 
single system-wide scattering matrix; (|3.6I) determine the 
boundary conditions and apply them to the system-wide 
scattering matrix to calculate the distribution function ex- 
pansion coefficients; Q3.7P calculate the distribution func- 
tion values within the spin valve; and (|3.8p calculate the 
spin density (spin accumulation), spin current, and spin- 
transfer torque using the distribution function. 



/>, k) = f°a + f x a x + Pa y + f z a 
7 T (r,k) 



C/(r,k) 



/Kr,k) 



U^(r,k) (6) 3.1 Discretization of the Boltzmann equation 



where a x , <J y and a z are Pauli spin matrices, cto is a 2 x 2 
identity matrix, and U is a unitary rotation matrix. The 
rotation matrix allows spin direction to point in an arbi- 
trary direction over the Fermi surface and as a function of 
position. 

In strong ferromagnets, the rotation matrix in Eq. ((SJ) 
is independent of the position on the Fermi surface or k 
as in Eq. This constraint is a consequence of the large 
difference in the Fermi surfaces for majority and minor- 
ity electrons in strong ferromagnets. The constraint arises 
because any transverse spin accumulation in the electrons 
near the Fermi surface rapidly dissipates. The transverse 
spins precess in the large exchange field and do so at differ- 
ent rates because of the complicated Fermi surfaces. The 
precessing transverse components rapidly dephase with re- 
spect to each other, leaving no net transverse moment. 



3 Formal technique 

Figure [1] shows a schematic picture of a spin valve - 
a magnetic multilayer structure where a thin-film non- 
magnet (spacer) is sandwiched between two thin-film fer- 
romagnet (FM) layers. The latter connect to reservoirs 
with non-magnetic leads. The main purpose of this paper 
is to solve the Boltzmann equation numerically for the dis- 
tribution function in such a spin valve structure. In reality, 
the cross sectional dimensions of these structures are much 
larger than typical layer thicknesses so that the transport 
in the interior of the sample is more important than that 
near the edges. The reservoirs have much larger cross sec- 
tional areas and hence smaller resistances than the active 
structures so that the transport is largely perpendicular to 
the layers. This combination of features suggests a simple 
model in which the layers are treated as translationally 
invariant in plane so that properties only vary perpendic- 
ular to the planes. Electrons move in all three directions, 
but the net current and all variation is in the cc-direction, 
i.e. /(r,k)=/>,k). 

The calculation proceeds in eight steps: (|3.1[) discretize 
the spin dependent Boltzmann equation using a numeri- 
cal mesh of the Fermi surface; (13. 2(1 solve the discretized 
Boltzmann equation for the eigensolutions in the non- 
magnetic and ferromagnetic bulk; (|3.3p use the eigenso- 
lutions to construct the layer scattering matrix for each 



To discretize the Boltzmann equation, we need a numer- 
ical mesh for the electron wave-vector k that can accu- 
rately and simultaneously describe all Fermi surfaces for 
both spin types and for all layers. A simple method for 
choosing a mesh is as follows. Choose a uniform mesh in 
the direction parallel to the interfaces (perpendicular to 
x): kjj. For each material, there could be several points 

on the Fermi surface that have the same kjj. We label 
their longitudinal wave- vectors fc™. A complete mesh for 
the Fermi surface is k^ = (A™, kjj). The mesh weights are 
determined by the area on each Fermi surface associated 
with fc™. This mesh may converge slowly, so it may be nec- 
essary to refine the kjj mesh to include a higher density of 
points. Let us assume one such mesh has sampling points 
{ki}^ with weighting factor {wi}^ . Using this mesh, the 
integration for any continuous function /i(k) on the Fermi 
surface can be discretized as: 



N 

/ fr(k)dk = y~V/i(k 



(J) 



We discretize the integrations in Eq. (g]) using this 
mesh. Assuming the system is one dimensional in x direc- 
tion, we have a discretized form of the Boltzmann equa- 
tion: 

" /; eE x = Y} V ~^W '// ', (8) 



dx 



where the subscript i and j mean that k is evaluated at 
ki or kj, for instance f? — f a {ki). The velocity matrix 
V and the scattering matrix B are 2N x 27V matrices (N 
from ki index, 2 from spin index) with matrix elements 



'j 



rtcrcr 



= v ax S aa 

* V 



WjP° j 6 atT > 



(9) 



5f? 



+ w i pfAi 



6 aa 



where 5?? — 5ijS a(T ' , 1/rf = J2j w jPiji IS the spin-dependent 
non-spin-ffip scattering rate, and 1/rf = J2j w jPij is the 
spin flip scattering rate. See Eq. The matrix B is sin- 
gular because ^2 ja , B? a =0. 
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3.2 Solution of the Boltzmann equation in the bulk 

In this section, we briefly summarize the results of Rcf . 26J 
and then generalize them to treat systems with non-collinear 
magnetizations. According to Ref. [26], the Boltzmann 
equation Eq. (JSj) can be solved by breaking it into a par- 
ticular equation 



E 



[V-^B)lfff = -eE x , 
and a homogeneous equation, 

j,cr' 

The solution to the particular equation Eq. (|10[) is 
FZ{xM) = -eE x Y}.B- 1 V]lf. 



(10) 



(11) 



(12) 



Due to the singularity of B, B X V is not the inverse of 
V- X B. Instead, the inverse matrix is defined as in Section 
2.9 of Ref. [27] using a singular value decomposition. Since 
the constant vector is perpendicular to the null space, the 
solution is well defined. Physically, the solution is the dis- 
tribution function of the current in bulk material. Since 
this distribution is well defined physically, we expect it 
to be mathematically as well, if we have formulated the 
problem correctly. 

Most solutions to the homogeneous equation Eq. (jlip 
vary exponentially in space: 

F„ CT (ar,k i )=^(k i )e A » a: with n€ [3,27V], (13) 

where g^(ki) = and X n are the n-th eigenvector and 
eigenvalue of V^B 



E 

3,"' 



3<T 



(14) 



See the Appendix in Ref. ^26] for more details. Since the 
bulk is translationally invariant in x direction, half of the 
eigenvalues A„ are positive and half are negative. The ma- 
trix V~ 1 B is defective, which means that the degenerate 
zero eigenvalue only has one eigenvector. Because it is de- 
fective, the homogeneous equation has two solutions that 
do not have exponential form 



F^(x,Vn) 



1 



xF[(x,ki) 



3,°' 



B 



(15) 



which can be verified by plugging them back in Eq. (|lip . 

Physically, F\ is the solution describing a uniform shift 
of the chemical potential and F2 describes a current car- 
rying solution having a spatially varying density and a 
associated uniform diffusion current. The rest of eigenso- 
lutions, Eq. (|T2"|) and Eq. (fT3]) . are exponential solutions 



that are necessary near interfaces. These solutions are not 
allowed in a uniform bulk because they diverge in one 
direction. Together, these solutions describe arbitrary so- 
lutions of the Boltzmann equation in each layer. 

Both the particular solution Fq and the homogeneous 
solution F2 are current carrying because of the k^-dependen 
of their summation terms. But, they describe current as- 
sociated with different processes. Fq describes the current 
due to the electric field, and F2 describes the current due 
to density gradients. Formally, F2 — A(E x x — Fq), where 
A is a constant, so Fo plus the electric field E x can be in- 
terchanged with F2 ■ Computationally, since we work with 
uniform current, we solve everything with E x — and 
no Fq and work with F^, Once we find the coefficient of 
F2, the physical solution is the corresponding Fq and E x . 
That is, we solve the problem as if the uniform current was 
exclusively due to diffusion with no electric field and then 
reinterpret the charge accumulation as an electric poten- 
tial and set the charge accumulation to zero. This inter- 
pretation makes sense due to the short screening length 
in metals. 

There is a natural quantization axis in a fcrromagnct 
defined by its magnetization. The distribution function is 
described by the eigensolutions in Eq. (fT5|) : 

/(a;,k i ) = /V T + /V i 

2N 



f° = Y,a°FZ(x,k i )=F° 



71=1 



(16) 



Here, ct t = (1/2)(I + a z ), a { = (1/2)(I - a z ), and the < 
are the expansion coefficients. 

In a non-magnet, there is no natural quantization axis, 
so we write the distribution function as in Eq. ([6]). We 
construct a different basis set F^, r = 0, x, y, z from F^ 
for the non-magnet. First, we know that the eigenvectors 
F£(x, ki) in Eq. (fT3|) and Eq. (fl~5]) break into separate 
eigenvectors for charge transport with F^ = F^ and spin 
transport with F^ = —F*. For instance, the eigenvectors 
with n = 1 and 2 correspond to charge transport because 
F± 2 = F\2' l n g enerai ; na H °f the eigenvectors (assume 
for the first half: n S [1, iV]) corresponds to the charge 
transport; the other half (for the second half: n e [N + 
1,2 AT]) is for the spin transport. Therefore, the basis set 
in a non-magnet can be constructed as: 



FT°{x,ki) = \ [Fj(.x,k0 +Fi(x,k 4 )] 



n G [1, 
(17a) 



fr^Or.k;) = - [Fl +N (x,ki) - F^ixM)} n e [1, 

(17b) 

Thus, in the non-magnet layers, the general solution for 
the distribution function is 

/>, k*) = fcjQ + fa x + Pa x + f z a z 



r 



N 

E 



a s n F^x,^)=F s - a s 



(18) 
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where s = 0,x,y,z, and a s n are expansion coefficients. 
Eq. (I17|) tells us that f x , f v , and f z share the same set of 
eigenvectors. This is because the Boltzmann solution does 
not depend on the choice of spin quantization axis in a 
non-magnet. 



3.3 Layer scattering matrix 

The eigensolutions are used to construct a scattering ma- 
trix which relates the boundary values of each layer. Fig- 
ure [2] shows a schematic picture of a spin valve, where 
lead/FM/spacer/FM/lead are labeled layer 1 to layer 5, 
and xq, 1,2,3,4,5 denotes the x coordinate of each interface. 
In Eq. (fl~6|) and Eq. (fl"8|) , we expand the distribution func- 
tion in each layer using a set of eigensolutions. From these 
expansions we construct the (layer) scattering matrix for 
layer m, § m , which relates the incoming electron distribu- 
tion functions f™ L/ , R and the outgoing ones f^L/n a * the 
left (L) and right (R) sides of the layer m (Figure [2]): 



f out 
J m,L 

f out 
im.R 



fin 
J m,L 

fin 
Jm,R 



(19) 



Here, the matrices in the ferromagnet (to — 2, 4) and 
non-magnet (to = 1, 3, 5) are respectively 



and fm'l = fl(Xm-l,K)> and fm,R = fl( X m,K ), 

k.f is the electrons' wave- vector on the Fermi surface with 
the (+/-) superscript indicating whether the electorn is 
right (+) or left (-) going. The superscripts "in" and "out" 
denote electrons going into or out of the layer. The defi- 
nition in Eq. (|20[) generalizes in a straightforward manner 
for "R" and "out." 

From Eq. p6[) , we can express the distribution function 
as f^(x,ki) = F a (x,ki) ■ a s in the ferromagnetic layers 
(m = 2, 4), then 



j.in,T An, I 
J m,L ' J m,L 

fin,0 fin,x fin,y fin,z 
■/ m,L ' Jm,L' ■> m,L > Jm,L 



(20) 



etc. 



j*out 
J m,L 

£OUt 

Jm,R 



fin 
J m,L 

fin 
im,R 



,K) 



F T (ir 

FJ-(a; m _i,k; 
FT(x m ,k+ 

P T (a; m _i,k; 
FJ-(a; m _i,k; 



(21) 



where and F™* are both 2N x 2N square matrices. 
Therefore, Eq. (fT9|) implies that the layer scattering ma- 
trix for ferromagnetic layers (m = 2, 4) is 



§m. — F° 



[f; 



mj 



(22) 



The layer scattering matrix for non-magnetic layers (m = 
1, 3, 5) is constructed in a similar way using Eq. fTS]) as 
the expansion and Eq. ([20| for m — 1, 3, 5. 



3.4 Interface scattering matrix 

Using the layer scattering matrix, we are able to relate 
the distribution function values at the two sides of a bulk 
layer. Next, we find an interface scattering matrix which 
connects the distribution functions across an interface. 
Right at the interface, the Boltzmann equation is not valid 
because the material properties vary rapidly. The distri- 
bution functions across the interface are related through 
the scattering matrix for the electron wave- functions at 
the interface. 

Like the layer scattering matrix in Eq. (|19[) . the in- 
terface scattering matrix R relates the incoming to the 
outgoing distribution functions. If the interface is specu- 
lar, when electrons scatter from the interface, the compo- 
nent of the wave vector that is parallel to the interface 
is conserved. In this case, the interface scattering matrix 
is block diagonal, with non-zero elements only for those 
states with the same parallel wave vector. On the other 
hand, defect scattering at the interface couples states with 
different parallel wave vectors and the interface scattering 
matrix becomes dense. In the present work, we neglect de- 
fect scattering at the interfaces for simplicity. Interfacial 
defect scattering has been considered by several authors 
[28129 30J. Xia et al. [30] treated interdiffusion at Co/Cu 
interfaces with first principles calculations and found that 
such defects caused only minor changes in the average 
transport properties for these interfaces. 

Consider an isolated interface between a non-magnet 
and a ferromagnet (NM/FM) and choose the spin quan- 
tization axis to be parallel to the magnetization of the 
ferromagnet. The interface between Lead 1 and FM2 in 
Figure [2] is such an interface if those two materials were 
extended to infinity. Suppose the wave-function for an 
electron on the Fermi surface in the corresponding NM is 
<fi(ki), which is orthonormal to other states on the Fermi 
surface: (0(ki)|0(kj)) = Sij. The wave-function for an 
electron on the Fermi surface in the FM is Vv(ki), which 
is also orthonormal to other states on the Fermi surface: 
(Vv(k;)|?/v(k,-)} = SijS, 



In the non-magnet, the wave-function for an electron 
moving toward the interface with spin pointing in an arbi- 
trary direction is written as a linear combination of spin- 
up and spin-down components: 



a<Mk+) 
b<P(K) 



(23) 



where a and b are the coefficients of the up and down 
spinor components. This incident state is scattered at the 
NM/FM interface, and the scattered states are 



■Pre 



E 

3 

E 

3 



K Na 0(k7) 
< f 'Vt0s + 

< F -Vl(k+ 



for x < 0, 



for x > 0, (24) 



where R^' 17 and T^ F:<T are the reflection and transmis- 
sion amplitudes for electron from k,; to kj for spin-up and 
spin-down electrons: a =t, I- 
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Fig. 2. Scattering matrices. The center of the figure lists the incoming and outgoing distribution functions at each of the layers. 
The top lists the layer, S n , and interface, R„ scattering matrices that relate the distribution functions. 



The 2x2 matrix distribution function is then denned 
by the outer product of the spinor coefficients. For in- 
stance, for the incident state 



fin — 



aa* ab* 
ba* bb* 



(25) 



Straightforward algebra reveals 



/ref(kj ) — R^ N / in (k+)R^ N 



/tr(k+) 



T^7i„(k+)T 



NF 

ij ' 



(26) 



where the reflection matrix (NM to NM) and transmission 
matrix (NM to FM) are 



XN 



R 



rpNF 

ij 



i?,T< T 



R 



NN,| 

"ij 



T NF, T q 

n T NF,1 
ij 



(27) 



Considering the electrons incident onto the interface 
from both sides of the interface, the scattering relationship 
Eq. (|26|) becomes 



/ N (k7) = J2 IC VN(k+)R^ + J2 T™ T /F(kr)T 



FNt i_/i,-\r rl FN 

ij ' 



(28a) 



/ F (kt) = £ R^ Ft / F (k-)R^ F + £ T^ Ft /N(k+)X 



NF 
ij " 



(28b) 



The matrix forms of the distribution functions / in 
Eq. (|2"5|) are expanded using cr-^ in the FM or Pauli ma- 
trices <jQ :X ,y,z in the NM as in Eq. ([6]). We represent the 
distribution functions by their expansion components fr 1 
in the FM and fo, x , v ,z in the NM as in Eq. After this 
transformation, the scattering formula Eq. (|28p is written 



fin 
JN 
fin 
JF 







r/°(x ,k+)i 


fn(xo,K ) 




/N(^o,k+) 


fU x o,K ) 


= R 


/^(x ,kt) 


/ N 0o,k 4 ) 




f§(x ,kf) 






fl{x ,kr) 






Jf( x o,K). 



f OLlt 

J N 
f out 



.(29) 

where R is the interface scattering matrix for the distri- 
bution functions across the interface, and the matrix ele- 
ments in R are obtained from the scattering matrices R 
and T in Eq. ((27)) . Note the apparent reversal of "in" and 
"out" in Eq. (gU) compared to that in in Eq. (H]). The di- 
rections in and out are defined with respect to the layers, 
rather than the interfaces and the incoming distribution 
for one of the layers is the outgoing distribution for one of 
the interfaces. 



3.5 System scattering matrix 

With all the layer and interface scattering matrices §1,2,3,4,5 
and Ri,2,3,4 (see Figure [5]), we are ready to construct a 
system- wide scattering matrix T that relates the incoming 
and outgoing distribution functions near the left and right 
reservoirs. The scattering matrix T is obtained by joining 
all the layer scattering matrices and interface scattering 
matrices in order: T = Si — Ri — §2 — R2 — §3 — R3 — §4 — 
R 4 - § 5 . 

The joining procedure is as follows. First we join Si 
to Ri (see Figure [2]), where Si is the layer scattering ma- 



trix that covers the interval [a 



> • 



(layer 1), and 



the interface scattering matrix that covers [x 1 , x~\ 
interface at x\): 



*i is 
(the 



f out " 
Jl,L 
f out 
.il,R. 




fin 
</l,R 
fin 
.J2,L _ 





s^ L S^ R 
s* L s« b 





fin 
Jl,h 




fin 



f out 
Jl,R 
f out 
J2,L 



(30) 



Here we have subdivided the distribution vectors (dis- 
cretized functions) into two subvectors for the values on 
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1 

u 



T 



the left (L) and right (R) . The S and K matrices are corre- with the rotated scattering matrix 
spondingly subdivided into four submatrices labelled LL, 
LR, RL, and RR, where, for example, the LR submatrix 
connects incoming values on the right to outgoing values 
on the left. We denote the joint scattering matrix by Ty, 
where the subscript "1" denotes a lead, "i" denotes an 
NM/FM or FM/NM interface, "f" denotes a ferromagnet, 
and "n" denotes the non-magnetic spacer layer. Ti; covers 
[xq",x^~] and relates /™/ c 



1 

ut 



(36) 



lifi 



a/out j /.in/out 
and J2.L '■ 



is joined with S3 using the same procedure as in 
(j3"2"|) and becomes Tufi n . Continuing by joining T n fi n 
3 , §4 , K4 and S5 gives a system wide scattering ma- 
trix T = Tiifinifii that relates the distribution functions 
near the left reservoir and the right reservoir: 



Eq 
to 



out 
L 



ft 

fin 
J2.L 



fin 
Jl,L 

f out 
J2,L 



(31) 



out 
L 



f Ol 
Jl, 
f out' 

■/5,R 



T 



fin 

/l,L 
fin' 
/5.R 



T 



/i(x ,k+) 



By eliminating the intermediate distribution functions /Jp 
in Eq. ([30|) . we have 



qpRL qpRR 



(32) 



where 



in/out (37) 

Eq. ([37]) is a condition on the boundary values of the dis- 
tribution functions for an arbitrary solution of the Boltz- 
mann equation in the multilayer. 

In the case that the two leads (layer 1 and layer 5) are 
semi-infinite, we need a system scattering matrix T that 
covers the interval fx7,a;j~l: 



r|pLL 

ijiLR _ 
ijiRL _ 

T RR — 



S^ R [1 + R^ L (1 - § rr ]RLL)-i§RR] K lR 
M RL (1 — S RR R^ L ) _1 S RL 

RL/-I sRRujLL^-IsRRtoLR 



(33) 



The scattering matrices described here and the method 
of joining them is more complicated than approaches us- 
ing transfer matrices, which relate the boundary values 
from one side to the other rather than outgoing to incom- 
ing boundary values. However, transfer matrix approaches 
become unstable as layers become thick and the present 
method does not. 

Using the same procedure as above we can join Ti; with 

52 and then with R2 to have a scattering matrix Tujg which 
covers [j;J,x^]. We continue to construct T U f, Tii n , Tnfi n , 
Tiifini, Tii fin if, Tufinifi, and Ti inninl . When coming to the the 
spacer layer, the eigensolutions and the scattering matrix 

53 use z'-axis instead of the z-axis (which is used for Tnfi) 
as the spin quantization axis (see Figure[2]). Therefore, we 
have to make an axis rotation at the spacer layer when 
joining Ti; n with S3. 

The distribution functions that Tufi and S3 relate are: 



f out 
7l,L 
fin 
J3,L 

f out' 
■>3,L 
f out 
^3,R 



= T 



nii 



fin 
/l,L 

f out 
/3,L 



fin 
J3,L 

fin' 
J3.R 



(34) 



The primed distribution functions are written using the z'- 
axis as the spin quantization axis. To join Tu n with S3, we 



^in/out 

write f' 



r p in/out /-in/out 

m terms of / 3 £ : / 3 £ 



ut/. 



in/out 
3,L 



u 



is the component representation [Eq. (|20[) ] of the unitary 
rotation matrix U in Eq. ©: /' = UfW => /' = U/. 
Then 



f out 

Jl.h 



= T 



nii 



fin 
Jl,L 



f out 
J 1,1. 



'3,L 



fin 

f out' 
^3,L 



(35) 



= T 



(38) 



Eq. (J3TJ) and Eq. (J38J) each have a total of 8N un- 
known coefficients in the expansions of /1 and f$ but only 
47V equations. To solve this problem, we study the bound- 
ary conditions, which apply restrictions on the distribu- 
tion functions /1 and /s, and so reduce the number of 
unknowns in the expansions. 



3.6 Boundary conditions 

By examining the properties of the leads and reservoirs, 
we can restrict the form of the distribution functions /1 
and /s in the leads, such that the number of unknown co- 
efficients in the expansions equals the number of equations 
in Eq. (I37|) or Eq. (f38|) . Thus, we can uniquely determine 
the distribution functions from Eq. (|3"T|) and Eq. (IBH . We 
treat semi-infinite leads and finite leads differently, and 
separate the discussion into two cases: 1) the leads are 
semi- infinite: the left/right lead extends to the left/right 
infinitely, and 2) both leads have finite length before con- 
necting to reservoirs. 



3.6.1 Semi-infinite leads 

If the leads are semi-infinite, the distribution function f\ 
in the left lead (layer 1) includes only the exponential 
eigensolutions with A n > 0. Then from Eq. (fl~8|) . we have 

N 

A°(x,k 4 ) = alf?(x,ki) +F°(x,k i ) + J2 "M^k,), 

n— 3 
A„>0 



N 



(39a) 
(39b) 



n=l 

n + N>0 
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We set a® = lto fix the current, because F 2 is the current 
carrying term. Due to the x-translational invariance of the 
exponential eigensolutions, half of the eigenvalues A n are 
positive and half are negative for charge transport and 
spin transport, respectively. Therefore there are 2N un- 
known a% x ' v ' z in total, N/2 for each component of 0, x, y, 
and z. 

Similarly, the distribution function f' 5 in the right lead 
(layer 5) includes only the exponential eigensolutions with 
A„ < 0: 

AT 

ti°(x, ki) = ®F°(x, ki) + J2 ki), (40a) 

n— 3 

A„<0 

JV 

f n T {xM)= E PnKi^i)- (40b) 

n=l 
A„ +N <0 

where r = x,y,z. The choice (3® = sets the potential at 
the right lead to be zero. Eq. (140|) also has 2N unknowns, 
t = 0, x, y, z, in total. 

Plugging Eq. (J39J) and Eq. (j40]) into Eq. l[3"8j) . we have 
4iV equations with 4iV unknowns, the expansion coeffi- 
cients a„'s and /3„'s can be solved. This gives the distri- 
bution functions in the left and right leads. The distri- 
bution function values inside the spin valve are calculated 
from the boundary values using the appropriate scattering 
matrices, as will be shown later. 

3.6.2 Finite leads 

In actual spin valve samples, the leads are usually short, 
after which the sample connects to essentially bulk ma- 
terial, here referred to as reservoirs. Electrons that enter 
the reservoir are much less likely to scatter back into the 
sample than to stay in the reservoir until they are ther- 
malized. This means that the reservoirs behave as perfect 
absorbers. The distribution of the electrons coming out 
of the reservoir is characteristic of the bulk, independent 
of the distribution of electrons coming in. The connec- 
tion between the reservoirs and leads has been studied in 
more detail by Berger [31] and Hamerle et al. [32] ■ When 
the leads have finite length and are connected to electron 
reservoirs, the exponential eigensolutions have no singu- 
larities and all of them should be included. In such a case, 
the form of the distribution functions is constrained by 
the following properties of a reservoir: the electrons leav- 
ing the reservoir have a bulk-like distribution function and 
the electrons with arbitrary distribution function can be 
absorbed by the reservoir. Based on these two properties of 
reservoirs, we propose that the distribution function near 
the reservoirs should satisfy: (1) for the electrons going 
from the reservoir to the lead, the distribution function 
is bulk-like, namely /j n L and f^ n R are spin independent 
and have contributions only from F° and F® in Eq. ([TT|) ; 
(2) for the electrons going from the lead to the reservoir, 
the distribution function has whatever structure it wishes, 
namely and have contributions from all F® ,x,v ' z . 



s L s R 



Lead 1 


FM2 


SP 3 


FM4 


Lead 5 


/!°L 




1 3,R 




s out 
J 5,R 


f out 
J l.L 


f 


/3°R 




f in 
/ 5,R 






x' 





x x x x 2 x 3 x 4 x 5 x 

Fig. 3. Back-propagation matrices. 



Glf y ' z = F^' y ' z , (41a) 
G^*(:ro,k+) = 0, (41b) 
G°^(x 5 ,k-) = 0. (41c) 

Essentially, G p is obtained by using the linear combi- 
nation of Fi(xo,kf ) and F q (xo, kf) to cancel F p (xo,kf), 
and similarly G q is obtained by using F\ (x$ , k~ ) and F p ( x$ , k^~ ) 
to cancel F q (x$ , k~ ) . 

Gi t 2 and G p form the basis for the electrons that have 
bulk-like right-going behavior at x = xq, and G\^ and G q 
form the basis for electrons that have bulk-like left-going 
behavior at x = x§ . The distribution functions that satisfy 
the requirements (1) and (2) in the leads are constructed 
as the following: 

f°(x, ki) = a°G°(x, ki) + G° 2 (x, ki) + £ a° p G° p (x, k ( ), 

p 

(42a) 

f^ z (x, K) = £ a x ' v ' z G p ' v ' z (x, K), (42b) 
p 

fl(x,ki) =/3 2 °G°(x,k 4 )+E/? 9 °G°(x,k l ), (42c) 

i 

f%' v ' z (x, ki) = £ p*>V'*G^*(x, ki). (42d) 

9 

In these equations, = 1 fixes the current, and /3° = 
fixes the chemical potential at the right boundary to be 
zero. Since the indexes p and q each take (N — 2)/2 val- 
ues for the 0-component and N/2 values for the {x, y, z}- 
components, the total number of unknown coefficients in 
Eq. |g5J) is AN. Therefore, by plugging Eq. (g2J) into Eq. (13"7)) . 
the coefficients a*,/?£ can be determined. 



3.7 System distribution function 

Once we have the distribution function values at the bound- 
aries, either near the lead/FM interface for the semi-infinite 
lead case or near the reservoir for the finite lead case, we 



To determine the form of the distribution functions 
near the reservoirs, let us first use F®' x,y,z to denote the 
eigensolutions with positive eigenvalues: A p > 0; and use 
p0,x,y,z ^ Q denote the eigensolutions with negative eigen- 
values: X q < 0. We construct a new set of basis functions 
Q0,x,y,z usm g ii ncar combinations of the old basis func- 
tions F2< x *>* such that, 



f^ z (x ) k i ) = 
fl(x,ki) = 
f%' v ' z (x,ki) = 
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can calculate the distribution function values everywhere 
inside the spin valve using the scattering matrices. For 
example, assume the scattering matrices S L covers the in- 



terval 



x o j x 3 ] 



and § fl covers the interval 



fin 
Jl,L 

fin' 
J3,R 

f out' 

■>3,R 
fin' 
J5.R 



x 3 i X 5 



" rout " 
/l,L 
f out' 

.^3,R _ 






§lr" 

^RR. 




fin' 

fout' 
/5,R 




"s£, 


RR . 





(43) 



/3R U * can be solved from Eq. (|43| . but the equations 
in Eq. (|43[) are redundant, so we choose the half of the 
equations that use the incoming boundary values rather 
than the outgoing ones, i.e., 



f out' (j^i, fin 1 fin 

/3,R — s RLil,L "T s RR,/3,i?, 

fin' _ K-R f 
J3,R — S LL/3,R 



R fout' 1 55/?. fin' 
^LR/5,i?' 



(44a) 
(44b) 



From these equations, we can calculate Z™^ 011 ' • Similarly, 
we can calculate the distribution function value elsewhere 
using a different pair of scattering matrices § L and S R . 



3.8 Transport properties 

With the distribution functions in hand, it is straightfor- 
ward to calculate transport properties h(k) by integrating 
over the whole Fermi surface: 



h= I h(k)dk l . 
/PS 



(45) 



Using Eq. ([7]), the integrations for spin density and spin 
current are discretized as 



N 



spin density: n s m (x) = 'y]w i f^ l (x,k i ), (46a) 

i=l 
N 

spin current: j s m (x) = w^f^x, k t ), (46b) 



where s = 0,x,y,z for m = 1,3,5 (non-magnetic layers), 
and s =t, I for m = 2,4 (ferromagnetic layers). The spin 

frame: 



current at x — x 3 is written in the x'-z 



Q( x s)=jE( x 3 )x'+il(^3 



(47) 



where the jf (x 3 ) is the longitudinal piece parallel to the 
right FM layer's magnetization [z' direction), and j 3 {x 3 ) 
is the piece perpendicular to z'-axis. From Ref. [llj we 
know that the perpendicular spin current is absorbed at 
the NM/FM interface, therefore the spin-transfer torque 
acting on the right FM layer is 



N st =jf(aJi")*'- 



(48) 



4 Applications 



In this section, we apply our numerical method to calcu- 
late the spin-transfer torque acting on the right ferromag- 
netic layer of a model spin valve. We compare the results 
with equivalent calculations using the drift-diffusion ap- 
proach and Slonczewski's hybrid theory. First, we describe 
the approximations we make to simplify the calculation 
so as to focus on the differences between the different ap- 
proaches. 



4.1 Approximations and their rationale 



Compared with the drift-diffusion method and circuit the- 
ory, the most important feature of the Boltzmann method 
is its treatment of electrons in the bulk moving in dif- 
ferent directions. In circuit theory, the average over the 
Fermi surface, the spin accumulation, is used to charac- 
terize the electrons inside that node (here a layer) . In this 
treatment all of the electrons inside a node are effectively 
aligned. In the drift diffusion approximation, the distribu- 
tion is modeled by its first moment, the spin accumulation, 
and second moment with respect to velocity, the spin cur- 
rent. This allows for greater flexibility in describing the 
distribution function, but clearly if higher moments are 
important, neglecting them will lead to errors. In a Boltz- 
mann equation treatment of the transport the distribution 
function is allowed its full flexibility. 

To study the different approaches we consider a simple 
model in which we ignore the actual shape and/or size of 
the Fermi surfaces and assume that the Fermi surfaces in 
both the non-magnet and the ferromagnet (both spin-up 
and spin-down) are perfectly spherical and are the same 
size. We use different mean free paths to distinguish the 
differences between the electrons in the non-magnet and 
the ferromagnet: Zn for non-magnet, Zp for spin-up and 
spin-down electrons in ferromagnet. The choice of iden- 
tical and spherical Fermi surfaces makes finding a wave 
vector mesh particularly simple. Since the problem is az- 
imuthally symmetric, there are no contributions to the 
distribution function that are not azimuthally symmet- 
ric so that only polar variation need be considered. We 
choose Gauss-Legendre sampling for the polar direction 
with typically 40 mesh points. 

For the interface scattering coefficients in Eq. (|27|). we 
treat the case of ideal interfaces with no defect scatter- 
ing. Consider an electron with wave-vector k x incident on 
an interface. With the Fermi surface assumption made in 
the previous paragraph, the wave- vectors for the reflected 
and transmitted electrons are — k x and k x , respectively. 
The matrix elements of the reflection and transmission 
matrices R and T in Eq. (j27|) are calculated in Ref. [[33.] - 
To allow for finite and spin-dependent interface resistance 
in the equal-Fcrmi-surface model, we assume (^-function 
scattering at the interface to give the following transmis- 
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Table 1. Material parameters used in the Boltzmann calcula- 
tion. 



Parameter 


Material 


Value 


Units 


Reference 


I 


Cu 


110 


nm 


M 


kt 


Cu 


450 


nm 


ED 


/I 


Co 


16.25 


nm 


M 


l L 


Co 


6.01 


nm 


M 


ki 


Co 


59 


nm 


M 


a T 


Co/Cu 


0.051 




M 




Co/Cu 


0.393 




m 



sion and reflection probabilities: 



I4 N 'T = K'T 



(49a) 
(49b) 



where a =| or J., and fcf is a discretization of k x - The 
parameter a CT is proportional to the square root of the 
strength of the <5-function-like interface potential. This can 
be read off from the horizontal axis in Fig. 1 of Ref. [33] us- 
ing experimental spin dependent interface resistance data. 

We also use the relaxation-time approximation: P?j = 
P a = Aps/tv and P$ = P sf = A FS /T sf , where A FS is the 
area of the Fermi surface. In this limit, the current carry- 
ing eigensolution F% in Eq. (fl~5)) reduces to F%{x, k,*) = 
x + vfl^/vjr, where V is the mean free path for different 
spins and Vp is the Fermi velocity. 

Using the algorithm described above in Sec. [2] and the 
approximations discussed above (equal, spherical Fermi 
surfaces), we calculated the spin-transfer torque for the 
spin valve shown in Figure [TJ For the rest of this paper, 
we assume the non-magnetic and ferromagnetic layers are 
composed of Cu and Co, respectively. The results are quite 
similar if Cu and Co are replaced by other non-magnetic 
and ferromagnetic metal. The input values in the Boltz- 
mann calculation are listed in Table [TJ 



4.2 Results and Comparisons 

In this section we test the drift-diffusion approach and 
Slonczewski's hybrid theory by comparing the results with 
those found with the Boltzmann equation. The spin-transfer 
torque acting on the right ferromagnet layer (m layer in 
Figure [JJ in a spin valve can generally be written in the 
following form [lj: 

N s R t (0) = v(0)^ m x (rh x M), (50) 

where the cross product has magnitude of sin#. In Slon- 
czewski's hybrid theory, [15125] 



V(0) 



1- 



A + B cos 6 A- Boost 



(51) 



The parameters A, B, and q± are calculated using the ma- 
terial parameters and geometries as shown in Ref. [15J. An 
equivalent spin-transfer torque formula was obtained by 
Manschot, et al., [35] independently. 

As discussed above in Sec. 14. 1[ the drift diffusion ap- 
proximation assumes that the distribution function has a 
simple form consisting of a uniform expansion and a con- 
tribution proportional to the velocity. Thus, we can ex- 
pect that the drift-diffusion approximation breaks down 
when the variation of the distribution function over the 
Fermi surface is more complicated. There are three sit- 
uations where more complicated behavior is introduced. 
The first is when the transmission through the interface 
depends strongly on wave vector as is typically the case 
[37) . In the immediate vicinity of the interface, the wave- 
vector dependence of the transmission gives the distribu- 
tion function a complicated variation over the Fermi sur- 
face. This variation includes contributions that decay on 
the order of the mean free path, see Eq. (fl"3|) . If the inter- 
faces are separated by more than this length, the strong 
variation decays between the interfaces, and the two ap- 
proaches can be brought into agreement through an ap- 
propriate choice of an effective interface resistance in the 
drift diffusion approach. However, when the interfaces are 
closer, the interaction of these exponential contributions 
between interfaces complicate the transport. Evaluating 
the importance of these effects requires a calculation us- 
ing realistic band structures, which is beyond the present 
calculations. We instead evaluate the other two situations 
where such difficult calculations are not necessary. 

The second situation in which the distribution function 
has a complicated angular dependence is when the spacer 
layer is thin compared to its mean free path and the mag- 
netizations are not collinear. Figure [4] compares the an- 
gular dependence of the torque calculated with the drift 
diffusion approximation to the torques calculated with the 
Boltzmann equation. In these calculations, the reflection 
at the interfaces has been set to zero so that the compli- 
cations described in the previous paragraph do not play 
a role. Inset (b) in Figure 0] shows that the torques agree 
when the spacer layer is very thick, and inset (a) shows 
that when they are thin, there are significant differences. 
The main panel shows the variation of the maximum of 
the torque curves as a function of thickness. The difference 
between the curves gives the corrections due to the com- 
plicated angular dependence of the distribution function. 

The torques decrease with thickness for two reasons. 
The large length scale decay is set by the spin diffusion 
length. When the layer is thicker than its spin diffusion 
length, spin-flip scattering leads to a significant decrease 
in the polarization of the current that crosses from one 
side to the other. For spacer layers thinner than their spin 
diffusion length, but longer than their mean free paths, 
the polarization of the current depends on the ratio of the 
effective polarized resistance to the effective unpolarized 
resistance. For these structures, in which the ferromag- 
netic layers are thicker than their spin diffusion length, 
the polarization of the current decays roughly like one 
over the thickness of the spacer layer. 
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Fig. 4. Spin-transfer torque at the right interface of the spacer 
layer in a spin valve with semi-infinite leads. Solid curves are 
calculated from the Boltzmann equation, dashed curves are 
from the drift-diffusion method. The two insets show the an- 
gular dependence of the torque for two specific thicknesses, 
t$ = 10 nm (a) and 3000 nm (b). The main panel shows the 
thickness dependence of the maximum values as a function of 
angle for the two approaches. The legend gives the thicknesses 
of the layers (inf=infmite). 



The final situation, in which the drift-diffusion ap- 
proach is not adequate to describe the full angular de- 
pendence of the distribution function, is at the interfaces 
between the leads and the reservoirs. Typically, in the 
drift-diffusion approach, the spin accumulation is set to 
zero at this point and the spin current is allowed to vary. 
The argument is roughly that the large total density of 
states there compared to in the leads forces the spin ac- 
cumulation to be small. In Sec. 13.6.21 we described how 
the greater flexibility available in the Boltzmann equa- 
tion allows the implementation of boundary conditions 
that treat the reservoir as an absorber. In Figure [51 we 
show the differences that can result from the differences 
in the boundary conditions. Both calculations show that 
the angular variation in the torque depends strongly on 
the length of the leads. However, the Boltzmann equation 
results are not as sensitive as those from the drift dif- 
fusion calculation. In fact, the drift diffusion calculation 
gives both the parallel and antiparallel states as unstable 
for an asymmetric enough junction. 

The results described above show that the drift dif- 
fusion approach does not work when the layers are thin. 
Slonczewski [15125] , developed a simple hybrid theory that 
overcomes some of these difficulties. In particular, it treats 
the left going and right going electrons in the spacer layer 
separately. This overcomes the errors illustrated in Fig- 
ure [U The theory then treats the transport in the rest 
of the system with an approach closely related to circuit 
theory |21|22j . The result is an analytic expression for the 
torque, Eq. (|5 1[) . Here, we compare this hybrid theory with 
the Boltzmann equation to test its validity. In addition, 
we explore the systematic behavior of the spin-transfer 
torque as a function of the spin valve geometry. Figure [5] 
shows the angular dependence of the spin-transfer torque 
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Fig. 5. Spin-transfer torque at the right interface of the spacer 
layer in a spin valve. The left panel is calculated using the 
Boltzmann equation, the right panel using the drift-diffusion 
method. The legend gives the thicknesses of the layers. 
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Fig. 6. Spin-transfer torque at the right interface of 
the spacer layer in a spin valve with layer thicknesses 
5 nm/40 nm/ts/l nm/180 nm with ts = 1 nm, 80 nm, and 
160 nm. The solid curves are calculated from the Boltzmann 
equation. Solid circles are calculated by from the hybrid the- 
ory. The latter do not depend on £3. The inset shows the £3 
dependence of l/rj(9) for 9 = 0° and 180° for this geometry. 
The legend gives the thicknesses of the layers. 



acting on the right FM (Co) layer for a spin valve with 
geometry: 

Cu(5 nm)/Co(40 nm)/Cu(t 3 )/Co(l nm)/Cu(180 nm). 

The spacer layer thickness £3 varies from 1 nm to 160 
nm. The magnitude of the spin-transfer torque reduces as 
spacer layer thickness £3 increases. Features of the torque 
are discussed in Ref. [T5] , 

In Slonczewski's hybrid theory p. 5 25], scattering in 
the spacer layer is ignored. This means the spacer layer is 
treated as a thin film. In the case <3 = 1 nm, the spacer 
layer thickness satisfies the condition of the hybrid theory. 
If we fit the spin-transfer torque curve calculated from the 
Boltzmann equation using the spin-transfer torque for- 
mula Eq. ([50]) from the hybrid theory (see Figure [6] for 
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Fig. 7. Spin-transfer torque at the right interface of 
the spacer layer in a spin valve with layer thicknesses 
200 nm/i2/l nm/1 nm/160 nm with t 2 = 10 nm and 160 nm. 
Solid curves are calculated from the Boltzmann equation, solid 
circles are from the hybrid theory. The inset shows the ti de- 
pendence of l/r]{9) for 6 — 0° and 180°. The legend gives the 
thicknesses of the layers. 
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Fig. 8. Spin-transfer torque at the right interface of 
the spacer layer in a spin valve with layer thicknesses 
5 nm/40 nm/1 nm/1 nm/ta with t& = 10 nm, 80 nm, and 160 
nm. All solid curves are calculated from the Boltzmann equa- 
tion. The inset shows the ts dependence of l/r](d) for = 0° 
and 180° for this geometry. The legend gives the thicknesses 
of the layers. 



the fit), we find that the fitted interface resistance values 
agree with the experimental values within 15 %. This is 
very good agreement considering the experimental values 
themselves are accurate only within 10 % to 20 %. How- 
ever, if the spacer layer thickness becomes comparable to 
the mean free path in Cu, the torque curves (the solid 
curve in Figure [5] with — 80 nm and 160 nm) cannot 
be fit by the hybrid theory for any values of the interface 
resistances. 

The inset figure in Figure [H] shows how 1/77(0°) (solid 
line) and l/r;(180 o ) (dash line) vary with £3 in the Boltz- 
mann calculation. These quantities are related to the crit- 
ical current for initiating a magnetization switching: from 
parallel (P) to antiparallel (AP) Jp^ap oc 1/77(0°) and 
from antiparallel to parallel Jap^p °c 1/7/(180°). So the 
curves in the inset figure of Figure [6] also show that the 
critical currents vary almost linearly with the spacer layer 
thickness 7:3, and both curves have similar slopes. Experi- 
mental measurements show the critical currents increasing 
with spacer layer thickness [35]. 

We have seen in Figure [6] that Slonczewski's hybrid 
theory fails when the spacer layer is thick. The break- 
down of the hybrid theory is also seen in Figure where 
we show how the spin-transfer torque curve changes with 
the thickness of the left ferromagnetic layer £2- The in- 
put values in the hybrid theory here in Figure [Jj are the 
same as those used in Figure [U In the case 7j 2 = 10 nm, 
which is small compared to the spin flip length = 59 
nm in the ferromagnet, the hybrid theory and the Boltz- 
mann calculation agree with each other very well. When 
t2 = 160 nm, 7j 2 becomes comparable to or larger than l^ f , 
the hybrid theory starts to fail because an approximation 
of the hybrid theory does not hold when i 2 > 1^. This 
approximation assumes the spin currents at two sides of 
the thick ferromagnetic layer are equal: Q(xi) ~ Q{x2) 



(see FigureEJ. But in this case of ti > 1^, Q(x\) depends 
on t 2 in a non-trivial way. 

Next, we study a spin valve with geometry: 

Cu(5 nm)/Co(40 nm)/Cu(l nm)/Co(l nm)/Cu(t 5 ), 

where the right lead length t$ varies from 10 nm to 160 
nm. Figure [8] shows how the spin-transfer torque curve 
acting on the second (thin) Co layer changes when we 
vary £5. A second bump around 6 = 30° appears in Fig- 
ure [8j as ts becomes large. From the spin-transfer torque 
formula Eq. (|50[) and Eq. (|5ip in the hybrid theory, we 
see that the second bump corresponds to the g_ term in 
Eq. (|51[) . The value of q_ is typically close to zero and 
negligible, but it becomes prominent when the spin valve 
is highly asymmetric. By asymmetry, we mean that the 
left and right sides of the spacer layer have different spin 
dependent properties. For instance, for a spin valve with 
the geometry 

Cu(5 nm)/Co(40 nm)/Cu(l nm)/Co(l nm)/Cu(160 nm), 

the left side of the spacer layer has 5 nm Cu, and 40 nm 
Co, and two Cu/Co interfaces, which can be considered 
mostly ferromagnetic, because both Co and Cu/Co inter- 
faces have spin dependent resistances. However, on the 
right side of the spacer layer, there is only 1 nm of Co, 
while there are 160 nm Cu and two Cu/Co interfaces. So 
the 160 nm Cu dilutes the ferromagnetic character of the 
Co bulk and the Cu/Co interfaces and makes the right side 
of the spacer layer more like a non-magnet. This asymme- 
try of the spin valve - ferromagnet-like on the left and 
non-magnet-like on the right - leads to the emergence of 
the second bump in Figure [8j 
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5 Summary 

In summary, we developed a complete numerical algorithm 
to solve the Boltzmann equation in multilayer heterostruc- 
turcs using a scattering matrix formalism. This method 
solves the spin-dependent Boltzmann equation in a non- 
magnet and a ferromagnet and matches the bulk solutions 
using an interface scattering matrix for the distribution 
functions. The final solution for the distribution function 
is found by imposing boundary conditions, either from in- 
finite leads or from the electron reservoirs. Our interest 
in using this method is to calculate spin-transfer torque 
in a spin valve structure. The results were found to agree 
with the Slonczcwski's hybrid theory for geometries typi- 
cally encountered in experiments but not when layer thick- 
nesses become large compared to mean free paths. The 
drift-diffusion method agrees poorly with the Boltzmann 
calculation due to the extreme approximations it makes. 

One of us (J.X.) is grateful for support from the De- 
partment of Energy under Grant No. DE-FG02-04ER46170 
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